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ABSTRACT 

In the standard model for structure formation, bound objects originate from the 
gravitational collapse of small perturbations arising from quantum fluctuations with 
random phases. In other scenarios, based on defects, structures are seeded by localized 
energy density. In principle, it is possible to differentiate between these models on the 
basis of their statistical properties; only in the former case is the initial density field 
an almost-perfect random gaussian field. In this paper, we investigate the use of the 
trispectrum of the galaxy density field, which is the connected four-point function in 
Fourier space, as a discriminant between gaussian and non-gaussian models. It has the 
advantage of having only weak non-linear growth. We define a related statistic r which, 
as a test of the gaussian hypothesis, is independent of cosmology, the power spectrum 
and biasing, in real space, and which is, in principle, a measure of the departure from 
gaussian statistics. For galaxy redshift surveys, the statistic depends on cosmology and 
bias only through the potentially observable parameter (3. We compute the expected 
errors on the estimate of r, and demonstrate with numerical simulations that it can be 
a useful discriminant of models, with the important proviso that any bias is linear on 
large scales. Whether it is the most effective method is uncertain and depends on the 
nature of the departure from gaussianity. 

Subject headings: Cosmology: theory, large-scale structure of universe. Methods: ana- 
lytical 



1. Introduction 



Models of primordial fluctuations that generated the large-scale cosmological structures can 
be divided into two broad classes: gaussian and non-gaussian. The simplest versions of inflationary 



- 2 - 



cosmology predict almost perfectly gaussian initial fluctuations, but such fluctuations can arise 
under more general conditions from the central limit theorem. The popularity of gaussian models is 
also due to their mathematical simplicity and the possibility they present for analytical calculations. 
On the other hand, density fields generated by topological defects (Vilenkin 1985; Vachaspati 1986; 
Hill et al. 1989; Turok 1989; Albrecht & Stebbins 1992) have non-gaussian initial conditions; the 
same can be said about fields generated by some versions of inflation (Allen et al. 1987; Kofman 
& Pogosyan 1988; Salopek et al. 1989). Convincing evidence against gaussian initial conditions 
(GIC) would rule out many scenarios, and/or point us towards a physical theory for the origin 
of primordial fluctuations. Thus a test of the gaussian nature of the initial conditions is of great 
interest for cosmology. 

Microwave background anisotropies probe cosmic fluctuations at a time when their statistical 
distribution should be close to its primeval form. To date, due to the limited signal-to-noise ratio 
of existing data, no conclusive evidence about the gaussianity of the initial conditions has been 
reached (e.g., Heavens (1998); Ferreira et al. (1998); Hinshaw et al. (1994); Falk et al. (1993); Luo 
k Schramm (1993); Gangui et al. (1994); Luo (1994); Smoot et al. (1994); Kogut et al. (1996); 
Banday et al. (1999); Tegmark & Bromley (1999)). An alternative is to analyze the present day 
density field of the galaxy distribution. This approach is complicated by the fact that the density 
field we observe today has already undergone non-linear gravitational evolution, and the observed 
non-gaussian nature of the galaxy distribution on small scales may be entirely the result of non- 
linear gravitational clustering combined with biasing effects. This non- linear growth makes any 
intrinsic non-gaussian signal more difficult to detect, and leads one to the conclusion that generally 
study of the microwave background is likely to be more profitable in detecting non-gaussian features 
(Verde et al. 2000c). However, one cannot exclude the possibility that non-gaussian features arise 
on physical scales which are difficult to probe with the microwave background, or non-GIC might 
produce a gaussian Sachs- Wolfe effect (e.g., Scherrer & Shaffer (1995)), so large-scale structure 
studies may still have a role to play in the test of the GIC hypothesis. An alternative method 
is to study number densities of rare/high-redshift objects (e.g., Chiu, Ostriker & Strauss (1998); 
Robinson, Gawiser & Silk (2000); Willick (2000); Matarrese et al. (2000); Verde et al. (2000a,b)). 

There have been numerous studies to calculate the departures from gaussianity induced by 
gravity, and to set up a test of the gaussian nature of the initial conditions. There are essentially 
three different approaches: a) using N-body simulations, starting from GIC and several non-GIC, 
compare the resulting clustering properties of the evolved fields (e.g., Moscardini et al. (1991), 
Matarrese et al. (1991), Weinberg & Cole (1992)) b) investigate the topological properties of fields 
generated from gaussian and non-gaussian conditions and set up a comparison (e.g., Moscardini 
et al. (1990); Coles et al. (1993); Vogeley et al. (1994); Avelino (1997)) c) measure the moments of 
the density or the velocity fields and compare them with those predicted from a gaussian distribution 
or different non-gaussian models (e.g., Fry (1984); Fry & Scherrer (1994); Scherrer (1992); Catelan 
& Moscardini (1994); Catelan & Sherrer (1995), Gaztanaga & Fossalba (1998a,b), Kim Sz Strauss 
(1998); Bernardeau (1994)). 



-3- 



In this paper, we take the last approach, but work in Fourier space. We concentrate on 
the 4-point function, for the following reasons. Under rather general conditions, the two-point 
function carries no model-independent information about the gaussian nature of the field (Fan & 
Bardeen 1994). The bispectrum (3-point function) is more promising, as it is zero for a gaussian 
field, but generally non-zero. However, it grows at second-order in perturbation theory, and this 
signal will place limits on the accuracy with which one can identify a primordial component (Verde 
et al. 2000c). The trispectrum (4-point function) has the advantage that it grows linearly, with 
contributions from non-linear growth appearing only weakly in perturbation theory. This will be 
quantified in Sec. 2. One might hope, therefore, that the linear growth extends to relatively small 
scales. In addition, the possibility is open that non-gaussian models have no bispectrum, or even 
negative initial bispectrum (e.g., Moscardini et al. (1991)). 

One can also probe the 4-point function in real space (e.g., Luo & Schramm (1993); Lokas 
et al. (1995); Chodorowski & Bouchet (1996)), or look at a subset of configurations of the four-point 
function, through power correlations (Feldman et al. 1994; Stirling & Peacock 1996). An alternative 
way to approach mildly non-gaussian fields is based on the Edgeworth expansion (Amendola 1996). 

The main difficulty with these approaches appears if the smoothed galaxy distribution is related 
to the underlying mass distribution by a non-linear transformation. In this circumstance, it is 
probably very difficult to distinguish non-GIC from a non-linear bias, so one has to assume that, 
at least on large scales, the bias is linear and deterministic (cf., Taruya et al. (2000)). In any 
eventuality, there is a realistic possibility that current large galaxy surveys such as the Anglo- 
Australian 2-degree Field (hereafter 2dF; Colless (1996)) and the Sloan Digital Sky Survey (York 
et al. 2000) will place some constraints on any initial departures from gaussian behavior. 

This paper is organized as follows: In section 2 we review the relevant statistical properties 
of the density field through its n-point functions, and introduce the r statistic as a non-gaussian 
discriminant. Section 3 contains the practical implementation of the method on numerical sim- 
ulations. Redshift-space distortions are considered in Section 4 and conclusions are presented in 
Section 5. 

2. Statistics of the density field in the linear and weakly non-linear regime 

The statistical properties of the fractional overdensity field <5(x) = [p(x) — can be char- 
acterized by the n-point correlation functions or, in Fourier space, by the n-point spectra. If the 
fluctuation field is gaussian, the connected part of the n-point function vanishes for n > 3 (e.g., 
Bertschinger (1992)). Thus the two-point function, or alternatively the power spectrum, completely 
specifies a gaussian distribution. 

To linear order in perturbation theory an initially gaussian distribution remains gaussian; in 
particular the the connected n-point functions and n-point spectra of an initially gaussian dis- 
tribution are zero as long as linear perturbation theory holds. On the other hand, if the initial 
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conditions are non-gaussian, the n-point spectra in the linear regime are the primordial ones, scaled 
by the n* ft -power of the linear growth factor. Therefore it would be possible to detect primordial 
non-gaussianity with a measurement of non-zero connected linear n-spectra for n > 3. 

We concentrate here on the four-point function in Fourier space (the trispectrum), because, 
as we demonstrate below, second-order contributions vanish for gaussian primordial fluctuations, 
so we hope to be able to use linear theory to smaller scales. This is tested numerically. Moreover 
this quantity can be made independent of the cosmological model, the degree of linear bias, and 
the power spectrum at least on large scales as we will show below. 

The data which we will work with are products of four Fourier coefficients (in practice we use 
the real part of this) : 

D a = 5ki5k 2 5k 3 5k 4 - (1) 

The mean value of this is zero (by homogeneity), unless the four wavevectors form a quadrilateral. 
a labels the set of 4 wavevectors. 

Note that the trispectrum itself is strictly the Fourier counterpart of the connected part of the 
four-point correlation function only. For a zero-mean field with non-zero connected (subscript c) 
four-point function, 

(5 kl <5 k2 5 k3 (5 k4 ) = 

(<5 kl 5 k2 )(5 k3 5 k4 ) + (2 perms.) + ((5 kl 5 k2 <5 k3 5 k4 ) c (2) 

where 

(S ki S kj ) = (27r) 3 P(k i )5 D (k t + k j ) (3) 
(<5 kl <5 k2 <5 k3 <5 k4 ) c = (2vr) 3 T(k,)^(k 1 + k 2 + k 3 + k 4 ). 

Here P(k) is the power spectrum, T is the trispectrum, 5 D (k) is the three-dimensional Dirac delta 
function, the angle brackets indicate the ensemble average and (k«) is a shorthand for (ki, k 2 , k3, k 4 ). 

Note that (c> kl <5 k2 <5 k3 <5 k4 ) has a specific volume dependence when one works with a discrete 
rather than continuous Fourier transform. The gaussian part involves products of two Dirac delta 
functions, whereas the connected part has only one. For a volume-limited survey of volume V, 
the discrete Fourier transform changes these delta functions to 5 D — ► V/(2ir) 3 multiplied by a 
Kronecker delta, so any comparison between the size of non-gaussian and gaussian parts of D a is 
volume-dependent. Any statement about the relative importance of the terms must therefore be 
made with some caution. We have to stress here that this is a general feature of (<5 kl • • • <5 k „) that 
are evaluated in a finite volume: the relative importance of the connected and non-connected parts 
is volume dependent. 

In order to demonstrate the lack of contribution to {D a } from second-order perturbation theory, 
we expand the density field to second order as: 

<5(x) =5 (1) (x)+<5 (2) (x) (4) 



- 5 - 



where 5^ is the linearly-evolved overdensity. 5^ is 0(5^ 2 ) and, if initial conditions are gaussian, 
represents departures from gaussian behavior due to gravitational evolution (e.g., Goroff et al. 
(1986), Bouchet et al. (1992)). Strictly speaking the galaxy distribution is not a mildly non-linear 
field, but a highly non-linear field filtered on some smoothing scale. The operations of smoothing 
and evolution do not commute, but experiments (e.g., Matarrese et al. (1997), hereafter MVH97) 
show that perturbation theory works well if the smoothed field is not too non-linear. 

To second order in perturbation theory the Fourier counterpart of the connected four-point 
correlation function can be obtained by applying to eq. (2) the substitution 

* = 4? + ^ + ^ + ■■ (5) 

and retaining only terms oc e. Here e is simply a bookkeeping parameter that can then be set to 
unity and the expression for 5^ can be found e.g., in Fry (1984) and Catelan et al. (1995). Hence 



<4AA3<w = << ) 4 1 2 ) 4 1 3 ) + 

(4 1 1 ) 4 1 2 ) 4 1 3 ) 4 2 4 ) )+cyc.(4terms). (6) 

(2) 

There are 4 cyclical terms involving the second order 5^ whose expression involves products 
of two linear coefficients, so for an initially gaussian field, the second-order contributions to the 
trispectrum are products of five coefficients and hence vanish. The first non-vanishing contribution 
to the connected trispectrum due to gravitational instability is oc e 2 as is the next (third order) 
term in (5). In general, the trispectrum will depend on all orders of correlations present initially, 
but in the linear regime, the only contribution is from the initial trispectrum (T initia i). Specifically, 
if V(t) is the growth factor: 

T = Ti n itial V(t)\ (7) 

For GIC, the trispectrum is zero, but (D a )cic has a disconnected part 

(D a ) GIC = (2irfP(k 1 )P(k 3 )5 D (k 1 +k 2 )5 D (k 3 +k 4 )+ eye. (8) 

If primordial fluctuations are non-gaussian, then (D a ) ^ {DalGlC- F° r fc-vector configurations 
where all |k| are different, {D a )cic vanishes, leaving (D a ) = (2-7r) 3 T(ki, k2, k3, k.4)5 D . We can 
thus characterize the departure from gaussian behavior by 

(<5ki<5k2 < 5k 3 ( W 

= {2Tr) 6 P(k 1 )P(k 2 )5 D (k 1 + k 3 )5 D (k 2 + k 4 ) + eye. 

+ T^JPik^P^Pik^Pik^S ^ + k 2 + k 3 + k 4 ), (9) 
where we have introduced the quantity 

r(ki,k 2 ,k 3 ,k4) (1Q) 
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r is independent of the volume, the cosmological model, linear bias and redshift (or alternatively the 
power spectrum amplitude), but in general depends on the kj (i.e. the scale and the configuration 
of the quadrilateral formed by the 4 A;- vectors) . Here for simplicity we will use equilateral k- vector 
configurations; with this choice isotropy demands that r depends only on |k|, for a given shape. 
We can then use a likelihood method to estimate r(k), and deduce whether it is consistent with 
the gaussian result (r = 0). 



2.1. Likelihood 

Our treatment here follows that of MVH97. For simplicity, we would like to be able to adopt a 
gaussian likelihood (italics will refer to the form of the likelihood function, not the initial statistics) , 
which will be the a posteriori probability distribution for r if we assume a uniform prior: 

1 



£(r) 



(27r)-(detC) 



exp 



a/3 



(11) 



Here M is the number of data, which have means /x a (r) = (D a ), and that can be evaluated 
using equation 9. The covariance matrix is C a p = ((D a — [i a )(Dp — up)). If r = the likelihood 
analysis is effectively a x 2 analysis. 

In the general (i.e. non-gaussian) case we cannot justify the applicability of the gaussian 
likelihood mathematically from the central limit theorem, but numerical experiments (section 3) 
support its use when many modes are used. 

The covariance matrix involves the eight-point correlation function, computable from Wick's 
theorem and, for continuous and discrete fields by the methods of MVH97: 



(6l #8) total ~~ 




(6462) (6364) (6 5 6 e )(6 7 6 8 ) + .. 


. 105 terms 


+{Si 6 2 ) (S 3 5 4 5 5 ) (6 6 6 7 6 8 ) + .. 


. 280 terms 


+ (6 1 62)(6 3 6 4 6 5 6e)(6 7 6 8 ) + 


. 210 terms 


+(61626364) (6 5 6 6 6 7 6 8 ) + ... 


35 terms 


+(6\6 2 ) (6 3 646 5 6 6 6 7 6 8 ) + ... 


28 terms 


+(6i6 2 6 3 ) (646 5 6 6 6 7 6 8 ) + ... 


56 terms 


+(61... 6 8 ). 





(12) 



These expressions apply both to continuous fields and discrete point processes, and the averages on 
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the r.h.s. are all connected parts only. Under the hypothesis of GIC, as long as linear theory holds , 
all terms involving the connected m-point correlation with m greater or equal than two are non- 
zero only in the presence of shot noise. Moreover for equal-side configurations the terms involving 
(Si5 m S n ) are also zero, since one cannot make a triangle from any 3 of the wavevectors. For clarity 
we neglect shot noise in the main text, but the calculations have been performed including it. 

For GIC, the only non-zero terms come from the two-point function, i.e. the power spectrum. 
We can test the gaussian hypothesis by including only these gaussian terms in the covariance 
matrix, reducing the process to \ 2 analysis. This should also be a good approximation for mildly 
non-gaussian fields (cf. Heavens (1998) for the microwave background). With this covariance 
matrix, a measurement of r which is significantly non-zero would only rule out gaussian statistics, 
but would not necessarily recover the correct value of r, which would require use of the correct 
non-gaussian covariance matrix. 



2.2. A priori estimation of the error on r 

We can readily calculate the expected variance of our estimator for r. The calculation follows 
MVH97 so we only sketch the calculation here. Ignoring shot noise, for square configurations we 
have 

Ma = (D a ) = (2ir) 6 P 2 (5 D ) 2 + (2^frP 2 5 D (13) 

and the variance is: 

*a = (D a )=4P 4 [W 3 S D } 4 (14) 
since only 4 terms are non-zero out of the 105. Hence 

2 (PlnC ^ fiajr = 0)P 2 

a a 

The number of uncorrelated squares in a thin shell in /c-space of width 6 (ink) is irk?g5(h\ k) where 
g = V I (27r) 3 is the density of states. Considering contributions from all the shells to the continuum 
limit we obtain 

1 1 Fornax 

°; 2 = w4vL k3d{lnk) - (16) 

kmax is set by the breakdown of second-order perturbation theory and k m i n by the size of the 
sample. Therefore the error on r scales as y/V and the maximum signal-to-noise is obtained by 
splitting the volume into subunits. This procedure to reduce the error when dealing with higher- 
order spectral statistics is widely used in other fields such as signal processing 2 and fluid mechanics 



x As shown above, second order contributions vanish for gaussian primordial fluctuations, therefore linear theory 
breaks down when third-order contributions become important. 

2 Where it is referred to as segment averaging. 
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see e.g., Brillinger (1975), Brillinger & Rosenblatt (1967), Lii et al. (1976) and has been already 
applied in cosmology by Matarrese et al. (1997) and Verde et al. (1998). This appears to give 
something for nothing, but of course this is not true. The resolution of the paradox is given in 
section 2.3. 

The minimum size of the samples may be set by the breakdown of perturbation theory, or 
possibly by aliasing difficulties if there is significant power on scales larger than the box. 

Note that the analysis in this section, and therefore the scaling (16), holds only as long as the 
field is close to gaussian so that non-gaussian elements of the covariance matrix can be neglected. 
As we said before, if the field is close to gaussian we are justified in adopting a gaussian likelihood, 
however if the distribution appears to be highly non-gaussian this error calculation and volume 
dependence is no longer valid. Since effectively P(k) is estimated from the data themselves, in 
principle one should propagate any error on P(k) through into the error on r. However, since the 
cosmic variance for the trispectrum is much larger than that on the power spectrum, this effect 
should be negligible. We confirmed this with numerical experiments. 

The square configuration we considered above is not the only possible one. We could also 
consider other configurations in which all the /c-vectors in the quadruplet have the same modulus. 
These configurations can be parametrized by the angle 9 between the first two vectors and 90° > 
9 > 0°: in the square case 9 = 90°, in the 'degenerate' case 9 = 0°. The above error analysis holds 
for all 9 > 0°. For the 'degenerate' case of 9 = 0°, there are modifications to the number of cyclic 
terms contributing to the means and covariance matrix, and also to the number of independent 
quadrilaterals. In this case, we have fi a = 2P 2 [(2tt) 3 5 d ] 2 + (2tt) 3 tP 2 5 d , cr 2 = 24P 4 [(2ir) 3 5 D } 4 , 
and the number of uncorrelated quadruplets is 27rk 3 gd(lnk). These reduce the error by a factor 
^2 x 2 4 /24 = 1.15, if k max remains the same. 



2.3. Volume-splitting procedure 

Volume-splitting appears to improve the signal-to-noise. We show here in outline that it gives 
no more information than relaxing the fixed quadrilateral shape and including more modes. 

Consider the n point spectrum and imagine that we wish to compute the expected error on the 
amplitude A: (6^ ■■■&k n ) = AS5 D where S D oc V and S is some fixed function. The leading term 
in the covariance, neglecting shot noise, will be ~ P n (5 D ) n . As already seen before the density of 
states is g oc V. If we fix the configuration of the fc-vectors, then the expected error on A, for a 
thin shell in fc-space, will be given by: 

2 gS\5 D f 

A <X pn^Dy ( 17 ) 

which implies the curious result that a a oc y( n ~ 3 )/ 2 . However, if we allow the shape of the fc-vector 
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configuration to change we obtain: 



g (n-l ))g2((5 D V ; 



where the exponent n - 1 of the density of states comes from the presence of the Dirac delta 
function, i.e. the fact that the polygon has to close and therefore has n—1 degrees of freedom. In 
this case we obtain the (more intuitive result) that: a a oc 1/VV when all shapes are allowed for. 

In the case of the power spectrum (as in the example of Press et al. (1992)) the two possibilities 
coincide (n=l), for this reason splitting the volume does not improve the errors. Note also that 
in the bispectrum case (n=3) of Matarrese et al. (1997) and Verde et al. (1998), the error is 
independent of the size of the sample because the triangle shape is fixed. For the trispectrum 
(n=4) a a oc y/V '. Here, as in all cases, one has the choice between considering all shapes, or 
keeping the shape fixed and subdividing. 



3. N-body tests 

In this section, we perform tests with real-space results from N-body simulations. We analyze 
square configurations, and ensure that no wavevector occupies more than one square, so the covari- 
ance matrix is approximately diagonal (cf. MVH97 for the bispectrum). The simulation is a 128 3 
particle, 100 h^ 1 Mpc side box, CDM-like simulation from the Hydra consortium (Couchman et al. 
1995), with parameters Q = 1.0, A = 0.0, as = 0.64, T = 0.25, and GIC. Shot noise is completely 
negligible in all scales of interest (k < 1; here and hereafter k is in units of h/Mpc). 

The breakdown of linear theory for the trispectrum is not known in advance, but from MVH97 
one can expect that the leading-order corrections for the trispectrum (i.e. contributions from the 
third order in perturbation theory) should be small at least up to k = 0.55 h Mpc^ 1 . A likelihood 
computation (cf., section 2.1) shows that for this square configuration, linear theory for the 4-point 
function should be valid up to /c max = 0.67, which sets the upper limit for our further analysis, 
giving 1250 squares in total. This limit is also apparent in Fig. 1, which plots 3 (D a ). For this 
wavenumber limit, x 2 is shown as a function of r in Fig. 2: the minimum reduced \ 2 is 0.8 and 
r = (— 1.5 ± 5.5) x 10 4 is consistent with zero. The a priori estimation for the error of Eq.(16) 
yields a T = 5.6 x 10 4 . Note that the estimate of r is independent of linear biasing because of the 
presence of the same number of 8 factors in the numerator and denominator of Eq. 10. 

We have also analyzed another N-body simulation with very non-GIC. The initial conditions 



3 The notation here follows MVJ97: D a is a single measurement of a statistical quantity i.e. the trispectrum for 
a given (e.g., square) k configuration on a given scale (e.g., fc=0.5). (D a ) is the ensemble average of this quantity. 
Since we assume that the volume under consideration is a "fair sample" , then the average of all D a in the simulation 
on the same scale and with the same k configuration is an estimate of (D a ). 
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were set by applying the following mapping to a gaussian field: 

6 — ► 5 2 - (5 2 ) (19) 

(the so-called x 2 niodel for initial conditions). The linear power spectrum is a power law with 
spectral index n = — 1. The simulation we used has 80 3 particles, in a 200 h~ 1 Mpc side box, with 
eg = 0.5. In this simulation the shot noise is not negligible, but it can be included with the methods 
of MVH97. This family of initial conditions could arise from evolution through inflation of an 
isocurvature CDM model for structure formation and is motivated by the fact that it accommodates 
galaxy formation at high redshift (e.g., Peebles ( 1999a, b)). Despite the low normalization, in this 
case we find that for all fc-space shells with k < 0.67 the minimum \ 2 value for r is nonzero at high 
significance: considering all A;- vectors in the range 0.15 < k < 0.67 for equilateral configurations 
and in the range 0.15 < k < 0.75 for degenerate configurations, we obtain r eq = (9.2 ± 1.6) x 10 5 
and Tdeg = (6.4 ± 1.3) x 10 5 respectively. Thus we can reject the GIC hypothesis with confidence, 
but note that for this highly non-gaussian model the covariance matrix is not accurate, so we can 
not rely on the actual measurement of r being reliable. 

Incidentally, for this simulation, the bispectrum also differs in the mildly non-linear regime 
from the one that would have originated by gravitational instability from GIC. However, at the 
bispectrum level, there is degeneracy between biasing and initial non gaussianity: in practice it 
would not be possible to assess if there is substantial bias / anti-bias or if the initial conditions were 
truly non-gaussian. As we have mentioned, the trispectrum method is bias- independent, at least 
as long as the bias is not strongly non-linear. 



3.1. Subdividing the Volume 

We have tested the subdivision procedure with static simulations of a non-gaussian density 
field derived from a gaussian 5{x) 

S(x) — >5(x) + e5 2 (x) (20) 

where e is some parameter. For simplicity, we take a white noise power spectrum between k m i n and 
knmx- Details of the resulting 2-point and 4-point functions are given in the Appendix, but they 
are summarized as 

P(k) = P g (k) + 2e 2 P g (k) 2 V k /(2irf (21) 

where P g (k) is the power spectrum of the underlying gaussian 5 field, and Vk is a /c-space volume 
defined in the appendix. The 4-point function is 

<<Mk 2 <Mk 4 )conn. ~ (48e 2 P 3 3 + e 4 P^V k /(27rf)V. (22) 

We choose P g = 2(2vr) 3 in the range 0.1 < k < 0.837 and zero elsewhere, and e = 0.4. In this 
case, the non-gaussian terms in the covariance matrix are small in comparison with the gaussian 
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one and we ignore them. We split the volumes into cubes of side 20 and 24 h Mpc, and make 
repeated simulations to reduce errors. In all, we analyzed 460 of the smaller boxes and 266 of the 
larger boxes, so the total volume is the same in the two cases. The results are shown in Fig. 3. The 
expected value of r is 1040, and is recovered within the errors. Note also that the errors (160 and 
250 respectively) scale as expected oc ./V where N is the number of boxes, and the signal-to-noise 
is indeed higher for the ensemble of smaller boxes, as expected. 

4. Complications with a real survey 

In section 3 we showed that the trispectrum is a useful discriminant between GIC and non- 
GIC in a very idealized case where the field is unbiased, and the positions of the particles are 
known in real space. Note however that galaxy catalogues have an average number of galaxies per 
unit volume that is about two orders of magnitude smaller than the one of the simulation used, 
and use the redshift as a third spatial coordinate. The resulting redshift-space map of the galaxy 
distribution is therefore distorted and shot noise is significant. This is not a trivial issue since we 
are pushing into the mildly non-linear regime, which is significantly affected by non-perturbative 
contributions. 

In this section, we show how to deal with these effects. 

4.1. Redshift-space distortions 

It is a convenient approximation to split redshift-space distortions into two components (Kaiser 
1987): a large scale distortion responsible of the squashing known also as the 'bull's eye' effect and 
a small-scale radial smearing responsible for 'Fingers-of-God'. The large-scale effect on individual 
Fourier components of the density fluctuation is well described by multiplication by the Kaiser 
factor (1 + /3/x 2 ), where \i is the cosine of the angle between the /c-vector and the line of sight, and 
(3 = $7q 6 /6 with b the linear bias parameter. 

The small-scale effect is hard to treat exactly and could potentially erase the signal we are trying 
to detect; in fact virialized motions on small scales produce a radial smearing and the associated 
Finger-of-God effect contaminates the wavelengths we are interested in. A successful model that fits 
the power spectrum in numerical simulation reasonably well (e.g., Hatton & Cole (1998); Peacock 
& Dodds (1994)) assumes that the small-scale velocity field is uncorrelated with density and has 
an exponential velocity distribution. Although not rigorously theoretically-motivated, it has been 
shown in VHMM98 that this modeling of the small-scale velocity dispersion works well as an 
addition to perturbation theory. In the Fourier domain, and in the distant observer approximation, 
the exponential velocity dispersion gives a damping factor, that, combined with the boosting factor 
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of the large-scale effect (Kaiser 1987) gives (Peacock & Dodds 1994): 

^^ J Vi+fcw/2 ' (23) 

where a is the pairwise velocity dispersion of galaxies. With this model for redshift-space distor- 
tions: 

CI + Bu 2 ) A 

(D a )(k) (D a ){k,vi) = (D a )(k) (24 ) 

VrL=i,..,4U + fe^Vf/2) 

As in VHMM98 we allow <r to be scale-dependent (to fit the power spectrum - there is some 
observational evidence for this; see Hamilton & Tegmark (2000)) and we reject the fc-vectors aligned 
too closely with the line-of-sight (see Scoccimarro et al. (1999) for an alternative model). 

We have tested the model (24) for square configurations by performing a x 2 analysis for the 
parameter r on the unbiased redshift-space catalogue created from the N-body simulation with 
GIC. 

Following VHMM98 the covariance matrix is modified as follows. The power spectrum is 
first replaced by P(k) = P(k)(l + /3/i 2 ) 2 and the resulting covariance matrix is then divided by 

nj=i(l + k 2 o- 2 /J, 2 /2), where the index i runs over the 8 /c-vectors that form the two squares. 

The knowledge of the real-space power spectrum is required in this model: in a realistic appli- 
cation (like for example Sloan and 2dF) an accurate fit for the galaxy real-space power spectrum 
will be known. Thus for the present work we shall assume that the real-space power spectrum is 
known: only the velocity dispersion a needs to be determined, and this can be done by fitting the 
power spectrum. The limit of validity of the small-scale redshift-distortion model for the trispec- 
trum can also in principle be determined, although in this case the limiting factor is the breakdown 
of perturbation theory. A likelihood analysis of the redshift-space Fourier modes can give a in a 
similar manner to VHMM98. 

The result for the redshift-space analysis is shown in Fig. 4. Again the true value for r is 
well recovered within the errors: r = (-1.0 ± 5.5) x 10 4 (reduced X 2 value is 0.95) for the square 
configuration and r = (0.8 ± 4.5) x 10 4 , for the degenerate configuration in good agreement with 
the predicted error of 4.3 x 10 4 obtained as outlined in section 2.2 . 

As a test of the velocity dispersion, we consider the quadrupole-to-monopole ratio of (D a ), 
following its use with the power spectrum Hamilton (1992) and the bispectrum (VHMM98, Scoc- 
cimarro et al. (1998)). For the degenerate configuration all vectors have the same [i 2 , so the 
quadrupole-to-monopole ratio can be simply defined by Rt = {D a )^ /{D a )(°\ where the quadrupole 
and monopole moments of T are 



<A*> (0) = \j\D a ){k, l i)d l i. 



(25) 
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In Fig. 5 we show how the model (24) reproduces the observed Rt in the GIC simulation. We follow 
the same procedure as VHMM98, where a is allowed to be slightly scale-dependent, constrained 
to fit the power spectrum, and k 2 fj? < 0.3. It would be possible to define a similar quantity also 
for non-degenerate configurations, but in this case (D a ) in redshift space would depend on two 
different angles to the line-of-sight which leads to more complexity. 

For the non-gaussian isocurvature field of Section 3.1, we find a r which is inconsistent with 
zero: e.g., for degenerate configurations Td eg = (6.4±1.3) x 10 5 . For non-gaussian models, in general, 
the signal r will be shape-dependent, (as for example in the particular case illustrated in section 
5). This dependence on the /c-vector configuration might therefore hold some extra information or 
allow stronger detection of a non-gaussian signal. 

4.2. Expected performance from PSCz and Sloan 

The PSCz survey (Saunders et al. 1998) is the largest nearly-all-sky survey, containing around 
15000 galaxies. The space density of galaxies is not very high, so shot noise is important beyond 
about 40 h~ l Mpc for k = 0.67. Assuming the initial field is close to gaussian, we use the gaussian 
covariance matrix, but with shot noise included. 

Before being able to perform the trispectrum analysis on this survey however still there are a 
number of unresolved issues. The most important are the practical effects of the rapidly-varying 
selection function and the radial nature of the redshift-space distortions. In fact the technique 
described in section 4.1 assumes that the spherical nature of the distortion can be neglected and 
the sky can be considered flat (the so-called 'distant-observer approximation', or 'plane-parallel 
approximation'). We must emphasize that this approximation does not hold for the PSCz, and one 
really has to go beyond the plane-parallel approximation, by using spherical harmonics for example 
(e.g., Heavens & Taylor (1995)). 

In the case of the 2dF and Sloan surveys, due to the bigger volume available, the problem of 
the radial nature of the redshift-space distortions is much reduced, and the selection function issue 
can be tackled as illustrated in the appendix. The la error achievable on r is very encouraging: 
1 x 10 4 for Sloan, if we restrict the analysis to the linear regime (k < 0.3h) and divide the survey 
conservatively into 100 h^ 1 Mpc side cubes. This reduces to 4 x 10 3 if wavenumbers up to k = 0.7h 
Mpc -1 are included. The error can be reduced even further by considering smaller boxes. 

4-2.1. Meaning of non zero value for r 

Most inflationary models predict deviations from gaussianity of the form described by equation 
(26)-applied to the potential rather than the density fluctuation field; this kind of non-gaussianity 
affects mainly the bispectrum (e.g., Verde et al. (2000c)). An order of magnitude calculation shows 
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that for the non-gaussian models analyzed in Verde et al. (2000c), the minimum deviation from 
gaussianity that can be detected with the large-scale-structure bispectrum is at least a factor of 
a few smaller than the one detectable with the large-scale-structure trispectrum; in these cases 
an analysis of the cosmic microwave background bispectrum is thus more promising. However, as 
already mentioned, primordial non-gaussianity might arise on physical scales which are difficult to 
probe with the microwave background, the large-scale-structure bispectrum might deviate from the 
perturbation theory prediction due to the presence of bias or antibias, and non-gaussian models 
might have no primordial bispectrum, or even negative initial bispectrum. It is nevertheless clear 
that other plausible non-gaussian models such as the \ 2 model described in section 3 can easily be 
ruled out with a trispectrum analysis of Sloan/2dF data sets. 



5. Conclusions 

We have presented a method for using the trispectrum, the four-point function in Fourier space, 
to discriminate between GIC and non-GIC from linear large-scale structure data. The advantage 
of this method is that for mildly non-gaussian fields the analysis in real space is independent 
of the power spectrum normalization, linear biasing, and cosmology, and, as a spectral method, 
the covariance matrix is more simply computed, in contrast to n— point correlation functions. In 
redshift space, cosmology and bias enter only through the measurable quantity f3, we show how 
to deal with redshift-space distortions (in the distant-observer approximation), we include in the 
calculations effects of shot noise following the method presented in MVH97 (although for clarity the 
shot noise contribution is ignored in the text and reported in the appendix) and it is straightforward 
to model the effects of varying selection function using the method of MVH97 (see the Appendix). 
The non-linear evolution of the trispectrum is expected to be rather weak in perturbation theory, 
so our expectation is that linear theory should provide a good description up to scales where the 
field becomes significantly non-linear. This is born out (and quantified) in simulations. Fig. (1) 
illustrates this point. The equilateral configurations bispectrum for the same simulation considered 
here, agrees with second-order perturbation theory up to k ~ 0.55 (MVH97). Fig (1) shows 
linear perturbation theory is adequate for square trispectrum configurations up to k ~ 0.67. We 
parameterize the departures from gaussian statistics by introducing the related quantity r, which 
is the 4-point correlation function in units of sj ' P{k\)P{k2)P{k^)P{ki). This quantity, in specific 
cases, can give us a meaningful measure of 'non-gaussianity'. We compute the gaussian variance of 
r, which allows a test of the gaussian hypothesis. The error on r can be computed straightforwardly, 
and is found to be in good agreement with internal error from N-body simulations. For mildly non- 
gaussian fields, we can expect that the use of a gaussian covariance matrix is still quite adequate, 
and in such cases a measurement of r can reliably be made, subject to the requirement that the bias 
is linear on large scales. For highly non-gaussian fields, the gaussian hypothesis can be rejected, 
but the measurement of r will be unreliable as the covariance matrix we use may neglect important 
terms. 
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Application to data will probably have to wait for completion of the 2dF and Sloan surveys, 
as the distant-observer approximation is not a good one for the nearly all-sky PSCz survey. In 
addition, the high shot noise in the PSCz survey means that the expected errors on r are large 
enough to make it hard to distinguish between GIC and some significantly non-gaussian fields. The 
2dF and Sloan surveys should do much better. One could compute r from volumes of side 100 hr 1 
Mpc with a random error of less than one per cent for a x 2 model as the one considered here, and 
could therefore put tight constraints on non-gaussian models. It would be helpful to apply this 
method to defect models, but unfortunately none is available at the required stage of evolution. 

The major uncertainty in this analysis is the effect of bias. A non-linear bias term acts like 
a non-gaussian initial field, and, as we are using linear perturbation theory here, we cannot use 
polygons of different shape to lift the degeneracy, as done by MVH97 to lift the degeneracy between 
bias and gravitational evolution in the bispectrum. The best hope is to constrain a combination 
of initial non-gaussianity and quadratic bias, and to argue that if the measured value of r is zero, 
then it would require a conspiracy unless both effects were absent. In principle it is possible to 
distinguish these two effects by considering differences at higher order in perturbation theory, as 
an Eulerian non-linear bias is applied to the evolved field, whereas a primordial non-gaussian field 
of the same mathematical form is applied to the initial field, and thus the evolution is different. It 
is an open question as to whether any realistic survey would have the signal-to-noise to do this. 



APPENDIX 

We consider the transformation of a gaussian field as follows: 

S(x) — > 6(x) +e5 2 (x) (26) 

This can be seen as the first two terms of a Taylor expansion of any non-gaussian field originating 
as a local mapping from a underlying gaussian one. 

We choose a toy model for the power spectrum of the underlying gaussian field P g : a top hat 
function between some k minimum and maximum. The power spectrum for the resulting field will 
then be: 

P(k) = P g (k) + j^e 2 J P g (k')P g (\ k - k' \)d 3 k' ~ 

P g {k) + 2e 2 P g {k) 2 V k /(2iif (27) 

where V k is the volume where the integrand in the previous equation is non-zero. The relevant 
quantities for our analysis, substituting the Dirac delta function by V/(27t) 3 are the following: 

<<M k2 <Mk 4 >con„. ~ (48e 2 P 3 + e 4 P 9 4 Vfe/(2^) 3 )V, (28) 
SIGNAL = r = - — — — — ~ 

(<)kA 2 )(<!>k3<W 
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48e 2 P 3 + e 4 P 4 Vfc/(27r) 3 

(29) 



(P| + 4e2 P3y fc/(2vr) 3 + 4e^[V k /(2^W) 



NOISE ~ y/(27r) 3 y^ (30) 



1,3 _ p 
"'max n min 

The covariance matrix is dominated by the gaussian terms if (5^ 5k 2 <5k 3 ^k 4 )conn. *C (<5ki <5k 2 ) (<5k 3 #k 4 ) , 
but at the same time we want to maximize the signal-to-noise ratio. 

This considerations suggest that the best choice is: P(k) = 2(2ir) 3 if 0.1 < k < 0.837 and zero 
elsewhere and e = 0.4. This allows us to use the gaussian covariance matrix, but also to have a 
reasonable signal-to-noise for the non-gaussian component from a single box of size of a few tens of 
Mpc. 



Shot Noise 

In the application of the method to a real survey the shot noise might be important. Here 
we report how to modify the relevant expressions in the presence of shot noise. In what follows n 
denotes the average number of particles per unit volume, and the superscript d indicates that these 
results apply to discrete, point processes. We follow the generating functional approach of MVH97 
to find the following: 

(Wc^(2vr) 3 P{h) + \ S D 0si + kj) (31) 
(5id m S n ) d c — (27rf5 D (k l + k m + k n )x 

{Bl mn + I [P{h) + P(k m ) + P(k n )} + ^} , (32) 

where B[ mn denotes the Bispectrum, 

(5 5 p 5 q 5 r ) d c — ► (2vr) 3 5 D (k + k p + k q + k r ) x 

— [-B( 0+P ) gr + perm. (6 terms)] + 

[{Po+ P + q + eye. (4 terms) )+ 

Po+p + Po+q + Po+r] + 



n 3 



{6i...5 5 ) d c — ► (27r) 3 ,5 D (ki+k 2 + k3 + k4 + k 5 ) x 



(33) 
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[ 5 (i+2)(3+4)5 +perm. (15 terms)] + 
^2 [^(1+2+3)45 + Perm. (10 terms)] j + 



i [(-Pi+2 + eye. (10 terms)) + 



n 3 



+Pi + P 2 + P 3 + P 4 + P 5 ] + 
1 



n 4 



(<5i...<5 6 )f — (2vr) 3 ( 5 D (k 1 + ... + k 6 ) x 

[-Bi2(3+...+6) + P erm - ( 15 terms) + 
- B (i+2)(3+4)(5+6) + Perm. (15 terms.) + 
Pi(2+3)(4+5+6) + Perm. (60 terms)] + 
J_[p 1 + ...+p 6 + 

P\+2 + perm. (15 terms) + 

-P1+2+3 + perm. (10 terms)] + 
1 1 



n 5 



<$i...<S 8 >c — (2vr) 3 «5 D (k 1 + ...+k 8 ) x 

[Pi2(3+...+8) + Perm. (28 terms) + 

- B i(2+3)(4+...+8) + perm.(168 terms) + 
^i(2+3+4)(5+...+8) + perm.(280 terms) + 
-8(1+2) (3+4+5) (6+7+8) + perm. (280 terms) + 
5 (i+2)(3+4+5+6)(7+8) + perm.(210 terms)] + 

J_[P 1 + ... + P 8 + 

n 3 

Pi+2 + perm. (28 terms) + 
Pi_l_2+3 + perm. (56 terms) + 
Pi+2+3+4 + perm. (35 terms) + 



(34) 



(35) 



(36) 

Neglecting shot noise when it is not negligible has two effects. The main effect is to overestimate 
H. E.g., the generalization of (8) is T — > T$n, where 

(2TrfT SN = 
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(2^ 3 (P(k 1 )+^)(P(k 3 ) + k)* D Qti + k 2 )5 D (k 3 + k 4 ) 

+<<5i...<5 4 )^ (37) 

The second effect is that the errors are underestimated. In particular, even in the gaussian case, 
there are additional connected terms in the correlations arising from shot noise. 



Selection function 

The above analysis is valid for volume-limited samples, where the mean number density of 
galaxies is independent of position. In a realistic catalogue the presence of a varying selection 
function makes the mean density position-dependent and might even induce spurious detection of 
non-gaussianity. However it is possible to treat this effect accurately by the generating functional 
approach as in MVH97. If n(x) is the mean density of such a catalogue, following Feldman et al. 
(1994), we can define a working fluctuation field composed by subtracting from the real catalogue 
a synthetic catalogue with no clustering, but with the same selection function, and then weighting 
suitably the combination: 

F(x) = 7ts(x)[n(x) - cm s (x)] (38) 

where w(x) is an arbitrary weighting function, 7 is a normalization factor, n s (x) = n(x)/a, and a is 
the "dilution" of the synthetic catalogue. We will then consider the limit a — ► to avoid shot noise 
in the synthetic catalogue. The n-point correlation functions in Fourier space are easily calculated 
by considering F as the superposition of the process / = 7u>(x)n(x) and f s = -«7H)(x)n s (x). The 
generating functional for F will then be: Zp(J) = Zf(J)Zf s (—aJ). Ignoring shot noise, 

Zfsi-uJ) = -oa j d 3 xj s (39) 

where J s = 7^(x)n s (x) J2 m s m exp(— ik m • x), and the ansatz for the generating functional for the 
field / is: 

Z f [J] = 



cxp 



: / d 3 xj(x ) -\J d 3 xd 3 x'J(x)J(x')$ mi . (x, x') 



-| / d 3 xd 3 x'd 3 x''J(x)J(x')J(x'')^ n . (x, x', x") 
+±Jd 3 xd 3 x'dy'd 3 x l J(x)J{x')J{x') ■ ■ ■ J{x l )dln.(x, ■ ■ ;x l ) 

Following the same procedure as in Matarrese et al. (1997) it is possible to obtain the n-point 
function in Fourier space by differentiating the generating functional; shot noise can be easily 
included by modifying the generating functional . 

The r in presence of a spatially-varying selection function therefore becomes: 

r — A (41) 

J 22 
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where 



Iij = J tfW(xV(x) (42) 
In general the effect of a spatially-varying selection function can be summarized as follows: 



7 2 = VI22 (43) 
V _^ Im. . s D — > lNN (44) 

22 

(45) 



I-/ 2 ' (2^ 



1 / 



N(N-q) 



n q Inn 
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Fig. l.-The 4-point correlation function in Fourier space (disconnected and connected) 
{D a ) = (5kj . . . <5k 4 ) is approximated by averaging square configurations with binned wavevectors 
from a CDM-like N-body simulation (solid line). Also shown is the linear perturbation theory pre- 
diction for gaussian initial conditions i.e. fJ, a (r = 0) from equation (9), for a square configuration 
of wavevectors. Errors are errors in the mean for each bin. For this square configuration, linear 
perturbation theory breaks down around k = 0.67h Mpc -1 . 

Fig. 2.— Minimum x 2 analysis for the parameter r for the gaussian CDM-like N-body simu- 
lation. Only the 'square' configuration for the trispectrum has been considered here and 0.15 < 
k < 0.67 (in units of h Mpc -1 ). r is a measure of the connected part of the trispectrum; for 
gaussian initial conditions r = (vertical line). To produce this graph, we use the gaussian initial 
conditions covariance matrix to compute the likelihood (equation 11). See sections 2.1 and 3 for 
further details. 

Fig. 3.— The x 2 analysis for 460 volumes of 20/i -1 Mpc side (dot-dashed line) and 266 volumes 
of 24/1- 1 Mpc side (dashed line) for O.K k < 0.837 (in units of h Mpc" 1 ). The y axis is normalized 
so that l-<7 limits are at y = 1. The predicted value for r is r ~ 1040 (vertical lines). This illustrates 
the scaling with the volume of the error on r (see Section 3.1 for more details). 

Fig. 4.— Minimum \ 2 analysis for the parameter r from the redshift-space unbiased CDM-like 
N-body simulation (see text for details). Only the square configuration has been considered here. 
The value for the velocity dispersion parameter a ~ 980km/s although it is slightly scale-dependent. 
The range of k- vectors considered is k < 0.67h Mpc -1 . 

Fig. 5.- Quadrupole-to-monopole ratio for the trispectrum Rt for degenerate square config- 
uration of a GIC simulation. The continuous line is the theoretical Rt using the redshift-space 
distortions as in VHMM98 (see section 4.1 for details). 
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